##########################################################################################

library(data.table)
library(optparse)
library(ggplot2)
library(Seurat)
library(monocle)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 单细胞表达文件
    input_file <- "~/20231121_singleMuti/input/testis_combined.annotationCellType.Rdata"

    ## 高变基因的数量
    #variable_num <- 1000

    ## 输出
    out_path <- "~/20231121_singleMuti/results/monocole"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(input_file)

## 细胞顺序
cell_order <- c(5,4,8,9,17,10,7,1,2)

cell_level <- c("SSC","Differenting&Differented SPG",
    "Leptotene","Zygotene","Patchytene","Diplotene","Early stage of spermatids","Round&ElongateS.tids","Sperm"
    )

###########################################################################################

dat_diff <- subset( scrnat , seurat_clusters %in% cell_order )

## 提取前10%的高变基因对应的数量
variable_num <- round(length(rownames(dat_diff@assays$RNA)) * 0.1)

###########################################################################################
## 寻找高变基因

DefaultAssay(object = dat_diff) <- "RNA"
dat_diff <- FindVariableFeatures(object = dat_diff , nfeatures = variable_num)
ordering_genes <- VariableFeatures(object = dat_diff)

###########################################################################################
## 描绘高变基因推测的时序情况
dat_tumor <- dat_diff
cell_use <- rownames(dat_tumor@meta.data)
ALL_matrix <- as(as.matrix(dat_tumor@assays$RNA@data[,cell_use]),'sparseMatrix')
feature_ann <- data.frame(gene_id=rownames(ALL_matrix),gene_short_name=rownames(ALL_matrix))
rownames(feature_ann) <- rownames(ALL_matrix)
ALL_fd <- new("AnnotatedDataFrame", data = feature_ann)
sample_ann <- dat_tumor@meta.data
rownames(sample_ann)<-colnames(ALL_matrix)
ALL_pd <- new("AnnotatedDataFrame", data =sample_ann)
ALL_cds_all <- newCellDataSet(ALL_matrix,phenoData =ALL_pd,featureData =ALL_fd,expressionFamily=negbinomial.size())

## 估计size facotr，标准化细胞之间的mRNA的差异
ALL_cds_all <- estimateSizeFactors(ALL_cds_all)

## 离散度值可以帮助我们进行后续的差异分析
ALL_cds_all <- estimateDispersions(ALL_cds_all)

## 添加num_cells_expressed，可以用于过滤少数细胞表达的基因
ALL_cds_all <- detectGenes(ALL_cds_all)

## 将感兴趣的基因集合嵌入对象，后续操作都要依赖于这个list
ALL_cds_all <- setOrderingFilter(ALL_cds_all, ordering_genes)

## 降低数据的维度
# 慢
ALL_cds_all <- reduceDimension(ALL_cds_all, max_components = 2 , method = 'DDRTree')

## 根据发育轨迹对细胞进行排序 #reverse = T
# 慢
ALL_cds_all <- orderCells(ALL_cds_all)

out_file <- paste0( out_path , "/testis.monocle.Rdata" ) 
save( ALL_cds_all, file = out_file )

###########################################################################################
type <- "testis"

ALL_cds_all$cell_type <- factor(ALL_cds_all$cell_type , levels = cell_level)


## 黑色的点代码用于构造轨迹的差异基因；灰色是背景基因；
## 红色是根据计算的基因表达大小喝离散度分布的趋势
image_name <- paste0( out_path , "/plot_ordering_genes.",type,".pdf" )
plot2 <- plot_ordering_genes(ALL_cds_all)
ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)

## 拟时序的时间分布
image_name <- paste0( out_path , "/cell_type_trajectory",type,".pdf" )
plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="cell_type", cell_size=1)
ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)
image_name <- paste0( out_path , "/cell_type_trajectory",type,".divide.pdf" )
plot2 <- plot_cell_trajectory(ALL_cds_all, color_by = "cell_type") + facet_wrap('~cell_type', nrow = 1)
ggsave(file = image_name , plot = plot2,width = 15,height = 5)

image_name <- paste0( out_path , "/clusters_trajectory",type,".pdf" )
ALL_cds_all$seurat_clusters_character <- as.character(ALL_cds_all$seurat_clusters)
plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="seurat_clusters_character", cell_size=1)
ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)
image_name <- paste0( out_path , "/clusters_trajectory",type,".divide.pdf" )
ALL_cds_all$seurat_clusters_character <- as.character(ALL_cds_all$seurat_clusters)
plot2 <- plot_cell_trajectory(ALL_cds_all, color_by = "seurat_clusters_character") + facet_wrap('~seurat_clusters_character', nrow = 1)
ggsave(file = image_name , plot = plot2,width = 15,height = 5)

## 时序
image_name <- paste0( out_path , "/pseudotime_trajectory",type,".pdf" )
plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="Pseudotime", cell_size=1)
ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)


###########################################################################################
##  寻找拟时相关的基因（拟时差异基因）
ordergene <- ordering_genes 
Time_diff <- differentialGeneTest( ALL_cds_all[ordergene] , core = 1 ,
    fullModelFormulaStr = "~sm.ns(Pseudotime)"
    )

Time_diff <- Time_diff[,c(5,2:4,1,6,7)]
image_name <- paste0( out_path , "/pseudotime_ordergene.tsv" )
write.table(Time_diff , image_name , row.names = F , sep = "\t" , quote = F)



